Assignment #3

# SEE modeldata package for new datasets
library(tidyverse)         # for graphing and data cleaning
library(tidymodels)        # for modeling
library(stacks)            # for stacking models
library(naniar)            # for examining missing values (NAs)
library(lubridate)         # for date manipulation
library(moderndive)        # for King County housing data
library(DALEX)             # for model interpretation  
library(DALEXtra)          # for extension of DALEX
library(patchwork)         # for combining plots nicely
library(dbplyr)            # for SQL query "cheating" - part of tidyverse but needs to be loaded separately
library(mdsr)              # for accessing some databases - goes with Modern Data Science with R textbook
library(RMySQL)            # for accessing MySQL databases
library(RSQLite)           # for accessing SQLite databases
library(plotly)            # for interactive plots

#mapping
library(maps)              # for built-in maps
library(sf)                # for making maps using geom_sf
library(ggthemes)          # Lisa added - I like theme_map() for maps :)

#tidytext
library(tidytext)          # for text analysis, the tidy way!
library(textdata)          
library(reshape2)
library(wordcloud)         # for wordcloud
library(stopwords)

Put it on GitHub!

Here is my GitHub link.

Local Interpretable Machine Learning

We will use the King County housing data.

data("house_prices")

# Create log_price and drop price variable
house_prices <- house_prices %>% 
  mutate(log_price = log(price, base = 10)) %>% 
  # make all integers numeric ... fixes prediction problem
  mutate(across(where(is.integer), as.numeric)) %>% 
  select(-price)

Part 1

Choose 3 new observations and do the following for each observation: * Construct a break-down plot using the default ordering. Interpret the resulting graph. Which variables contribute most to each observation’s prediction? * Construct a SHAP graph and interpret it. Does it tell a similar story to the break-down plot? * Construct a LIME graph (follow my code carefully). How close is each original prediction to the prediction from the local model? Interpret the result. You can also try using fewer or more variables in the local model than Lisa used in the example.

set.seed(494)

house_split <- initial_split(house_prices, 
                             prop = .75)
house_train <- training(house_split)
house_test <- testing(house_split)

house_ranger_recipe <- 
  recipe(formula = log_price ~ ., 
         data = house_train) %>% 
  step_date(date, 
            features = "month") %>% 
  update_role(all_of(c("id",
                       "date")),
              new_role = "evaluative")

house_ranger_spec <- 
  rand_forest(mtry = 6, 
              min_n = 10, 
              trees = 200) %>% 
  set_mode("regression") %>% 
  set_engine("ranger")

house_ranger_workflow <- 
  workflow() %>% 
  add_recipe(house_ranger_recipe) %>% 
  add_model(house_ranger_spec) 

set.seed(494)
house_ranger_fit <- house_ranger_workflow %>% 
  fit(house_train)

house_rf_explain <- 
  explain_tidymodels(
    model = house_ranger_fit,
    data = house_train %>% select(-log_price), 
    y = house_train %>%  pull(log_price),
    label = "Random Forest"
  )

## Preparation of a new explainer is initiated
##   -> model label       :  Random Forest 
##   -> data              :  16210  rows  20  cols 
##   -> data              :  tibble converted into a data.frame 
##   -> target variable   :  16210  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.3 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  5.020731 , mean =  5.665231 , max =  6.696937  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.3328418 , mean =  0.0004361088 , max =  0.2467335  
##   A new explainer has been created! 

set.seed(494)
new_obs1 <- house_test %>% slice(4239)
new_obs2 <- house_test %>% slice(294)
new_obs3 <- house_test %>% slice(3)

10^(new_obs1$log_price)

## [1] 470000

10^(new_obs2$log_price)

## [1] 425000

10^(new_obs3$log_price)

## [1] 4e+05

Break-Down Plots

set.seed(494)

house_pp_rf1 <- predict_parts(explainer = house_rf_explain,
                       new_observation = new_obs1,
                       type = "break_down")

house_pp_rf2 <- predict_parts(explainer = house_rf_explain,
                       new_observation = new_obs2,
                       type = "break_down")

house_pp_rf3 <- predict_parts(explainer = house_rf_explain,
                       new_observation = new_obs3,
                       type = "break_down")

plot(house_pp_rf1)

plot(house_pp_rf2)

plot(house_pp_rf3)

• The average log price remains the same across all observations, as it is the average log price when applied to all of the training data. This log price is 5.665, or $462,381. We see that holding latitude constant at 47.698 affects the overall average price the most, and increases the log price by 0.084, or $98,666.95. The predicted log price for this observation is 5.653, or $449,779.90.

• The second plot shows us that fixing our latitude at 47.5197 decreases the average log price by 0.051, or $51,231.30. We also see that holding this latitude constant affects the price the most. Finally, the predicted log price for this observation is 5.599, or $397,191.50.

• The third plot shows us that fixing our latitude at 47.6127 increases the log price by 0.109, or $131,911.10. Holding this latitude constant affects the predicted price the most. This predicted log price is 5.59, or $389,045.10.

Shapley Additive Explanationss (SHAP) Plots

house_rf_shap1 <-predict_parts(explainer = house_rf_explain,
                        new_observation = new_obs1,
                        type = "shap",
                        B = 10
)

house_rf_shap2 <-predict_parts(explainer = house_rf_explain,
                        new_observation = new_obs2,
                        type = "shap",
                        B = 10
)

house_rf_shap3 <-predict_parts(explainer = house_rf_explain,
                        new_observation = new_obs3,
                        type = "shap",
                        B = 10
)

plot(house_rf_shap1)

plot(house_rf_shap2)

plot(house_rf_shap3)

• The first SHAP plot shows a similar effect to the first break-down plot, where grade is the highest contributing variable, however this time it is in a negative way. The second and third SHAP plots both show that latitude is the highest contributing variable, which is in accordance with the break-down plots, this time both affecting the predicted price in the same way as the break-down plots.

Local Interpretable Model-Agnostic Explanation (LIME) Plots

set.seed(494)

model_type.dalex_explainer <- DALEXtra::model_type.dalex_explainer
predict_model.dalex_explainer <- DALEXtra::predict_model.dalex_explainer

house_lime_rf1 <- predict_surrogate(explainer = house_rf_explain,
                                    new_observation = new_obs1 %>%
                                    select(-log_price), 
                                    n_features = 5,
                                    n_permutations = 1000,
                                    type = "lime")

house_lime_rf2 <- predict_surrogate(explainer = house_rf_explain,
                                    new_observation = new_obs2 %>%
                                    select(-log_price), 
                                    n_features = 5,
                                    n_permutations = 1000,
                                    type = "lime")

house_lime_rf3 <- predict_surrogate(explainer = house_rf_explain,
                                    new_observation = new_obs3 %>%
                                    select(-log_price), 
                                    n_features = 5,
                                    n_permutations = 1000,
                                    type = "lime")

house_lime_rf1 %>% 
  select(model_r2, model_prediction, prediction) %>% 
  distinct()

house_lime_rf2 %>% 
  select(model_r2, model_prediction, prediction) %>% 
  distinct()

house_lime_rf3 %>% 
  select(model_r2, model_prediction, prediction) %>% 
  distinct()

set.seed(494)

plot(house_lime_rf1) +
  labs(x = "Variable")

plot(house_lime_rf2) +
  labs(x = "Variable")

plot(house_lime_rf3) +
  labs(x = "Variable")

• Plot 1 shows us that the predicted log price is about 5.577, and that grade is once again the most important in the local model. This model also performs quite poorly, as the explanation fit is only 0.11. This is the exact prediction (to three digits) that the original model had.

• Plot 2 shows us that the predicted log price is about 5.86, and that the living room area and latitude are the most important in the local model, albeit in opposite ways. This model performs better than the first, as the explanation fit is 0.45. This is the exact prediction (to two digits) as the original model.

• Plot 3 shows us that the predicted log price is about 5.46, and that the grade is the most important to the local model. This model performs the worst out of the three, as the explanation fit is only 0.1. This is the exact prediction (to two digits) as the original model.

Part 2

Describe how you would use the interpretable machine learning tools we’ve learned (both local and global) in future machine learning projects? How does each of them help you?

• These methods provide a helpful alternative to a simple variable importance plot, as we are able to identify variable importance on the level of individual observations. We are able to see why each observation contains the result that it does. The break-down plot method allows us to see how the entire training dataset is affected when we apply a certain model to it. The SHAP plots allow us to see how each variable contributes if we change the order of the variables (in this assignment, the contribution changes because the random forest model is not additive). Through this, we are able to see how true the effect is through the boxplots on top of the SHAP plot. Finally, the LIME plots allow us to see variable importance in conjunction with model performance, which gives us more insight into how trustworthy our result is.

SQL

I will use the airlines data from the SQL database that Lisa used in the example in the tutorial.

Tasks

1. Create a SQL chunk and an equivalent R code chunk that does the following: for 2017, for each airport (with its name, not code), and month find the total number of departing flights, the average distance of the flight, and the proportion of flights that arrived more than 20 minutes late. In the R code chunk, write this out to a dataset.

con_air <- dbConnect_scidb("airlines")

DESCRIBE flights

Displaying records 1 - 10

Field	Type	Null	Key	Default	Extra
year	smallint(4)	YES	MUL	NA
month	smallint(2)	YES		NA
day	smallint(2)	YES		NA
dep_time	smallint(4)	YES		NA
sched_dep_time	smallint(4)	YES		NA
dep_delay	smallint(4)	YES		NA
arr_time	smallint(4)	YES		NA
sched_arr_time	smallint(4)	YES		NA
arr_delay	smallint(4)	YES		NA
carrier	varchar(2)	NO	MUL

DESCRIBE airports

9 records

Field	Type	Null	Key	Default	Extra
faa	varchar(3)	NO	PRI
name	varchar(255)	YES		NA
lat	decimal(10,7)	YES		NA
lon	decimal(10,7)	YES		NA
alt	int(11)	YES		NA
tz	smallint(4)	YES		NA
dst	char(1)	YES		NA
city	varchar(255)	YES		NA
country	varchar(255)	YES		NA

late_flights_over20 <- tbl(con_air, "flights") %>% 
  filter(year == 2017) %>%
  group_by(origin, month) %>% 
  summarize(tot_depart = n(),
            avg_dist = mean(distance),
            prop_late_over20 = mean(arr_delay > 20)) %>% 
  inner_join(tbl(con_air, "airports"),
             by = c("origin" = "faa")) %>% 
  select(name, month, tot_depart, avg_dist, prop_late_over20)

late_flights_over20_df <- collect(late_flights_over20)
late_flights_over20_df

late_flights_over20 %>% 
  show_query()

## <SQL>
## SELECT `origin`, `name`, `month`, `tot_depart`, `avg_dist`, `prop_late_over20`
## FROM (SELECT `origin`, `month`, `tot_depart`, `avg_dist`, `prop_late_over20`, `name`, `lat`, `lon`, `alt`, `tz`, `dst`, `city`, `country`
## FROM (SELECT `origin`, `month`, COUNT(*) AS `tot_depart`, AVG(`distance`) AS `avg_dist`, AVG(`arr_delay` > 20.0) AS `prop_late_over20`
## FROM `flights`
## WHERE (`year` = 2017.0)
## GROUP BY `origin`, `month`) `LHS`
## INNER JOIN `airports` AS `RHS`
## ON (`LHS`.`origin` = `RHS`.`faa`)
## ) `q01`

• With the dataset you wrote out, create a graph that helps illustrate the “worst” airports in terms of late arrivals.

late_plot <- ggplot(late_flights_over20_df %>% 
                    mutate(month = as.factor(month),
                    mean_prop = mean(prop_late_over20)) %>% 
                    group_by(mean_prop, name) %>% 
                    arrange(desc(mean_prop)) %>% 
                    head(200), 
                    aes(x = tot_depart, y = mean_prop)) +
geom_point(aes(size = avg_dist, color = name), alpha = .5) +
labs(x = "Average Distance Traveled",
     y = "Total Departures",
     title = "Which airports had the largest % of flights that \nwere more than 20 minutes late in 2017?") +
theme(legend.position = "none")

ggplotly(late_plot)

Create a table with 6 or fewer rows and 3 or fewer columns that summarizes which airport is the “worst” in terms of late arrivals.

late_flights_over20_simple <- tbl(con_air, "flights") %>% 
  filter(year == 2017) %>%
  group_by(origin) %>% 
  summarize(tot_depart = n(),
            prop_late_over20 = mean(arr_delay > 20)) %>% 
  inner_join(tbl(con_air, "airports"),
             by = c("origin" = "faa")) %>%
  ungroup() %>% 
  group_by(name) %>% 
  select(name, tot_depart, prop_late_over20, -origin) %>% 
  arrange(desc(prop_late_over20)) %>% 
  head(6)
  
late_flights_over20_simple_df <- collect(late_flights_over20_simple)
late_flights_over20_simple_df

2. An original SQL query and plot

Find the proportion of flights that arrive on time during winter months, with average flight distance, and sdee which airlines perform best by this metric.

on_time_winter <- tbl(con_air, "flights") %>% 
  filter(year == 2017, month %in% c(11, 12, 1, 2, 3)) %>%
  group_by(carrier) %>% 
  summarize(tot_depart = n(),
            avg_dist = mean(distance),
            prop_on_time = mean(arr_delay <= 0)) %>% 
  select(tot_depart, avg_dist, carrier, prop_on_time)

on_time_winter_df <- collect(on_time_winter)
on_time_winter_df

on_time_winter %>% 
  show_query()

## <SQL>
## SELECT `tot_depart`, `avg_dist`, `carrier`, `prop_on_time`
## FROM (SELECT `carrier`, COUNT(*) AS `tot_depart`, AVG(`distance`) AS `avg_dist`, AVG(`arr_delay` <= 0.0) AS `prop_on_time`
## FROM `flights`
## WHERE ((`year` = 2017.0) AND (`month` IN (11.0, 12.0, 1.0, 2.0, 3.0)))
## GROUP BY `carrier`) `q01`

on_time_winter_df %>% 
  ggplot(aes(x = prop_on_time, y = fct_reorder(carrier, prop_on_time), fill = avg_dist)) +
  geom_col() +
  labs(title = "Airlines with the Highest Proportion of \nOn-Time Flights during Winter Months",
       x = "Proportion of On-Time Arrivals",
       y = "Carrier",
       fill = "Average Flight Distance")

Function Friday

geom_sf() Tasks

states <- st_as_sf(maps::map("state", 
                             plot = FALSE, 
                             fill = TRUE))

counties <- st_as_sf(maps::map("county", 
                               plot = FALSE, 
                               fill = TRUE))

states <- states %>%
  mutate(area = as.numeric(st_area(states)))
head(states)

1. Change the color scheme of the map from the default blue (one option could be viridis).

ggplot(data = states) +
  geom_sf(aes(fill = area)) +
  scale_fill_viridis_c(option = "C") +
  coord_sf(xlim = c(-127, -63), 
           ylim = c(24, 51), 
           expand = FALSE) +
  theme_minimal()

2. Add a dot (or any symbol you want) to the centroid of each state.

ggplot(data = states) +
  geom_sf(aes(fill = area)) +
  scale_fill_viridis_c(option = "C") +
  stat_sf_coordinates(color = "white") +
  coord_sf(xlim = c(-127, -63), 
           ylim = c(24, 51), 
           expand = FALSE) +
  theme_minimal()

3. Add a layer onto the map with the counties.

ggplot() +
  geom_sf(data = states, aes(fill = area)) +
  geom_sf(data = counties, fill = NA, color = "black") +
  scale_fill_viridis_c(option = "C") +
  coord_sf(xlim = c(-127, -63), 
           ylim = c(24, 51), 
           expand = FALSE) +
  theme_minimal()

4. Change the coordinates of the map to zoom in on your favorite state.

ggplot() +
  geom_sf(data = states %>% filter(ID == "new york"), aes(fill = area)) +
  geom_sf(data = counties, fill = NA, color = "black") +
  scale_fill_viridis_c(option = "C") +
  coord_sf(xlim = c(-80, -71.8), 
           ylim = c(40.4, 45.1),
           expand = FALSE) +
  theme_minimal()

tidytext tasks

These are tweets from Twitter handles that are connected to the Internet Research Agency (IRA), a Russian “troll factory.” The majority of these tweets were posted from 2015-2017, but the datasets encompass tweets from February 2012 to May 2018.

Three of the main categories of troll tweet that we will be focusing on are Left Trolls, Right Trolls, and News Feed. Left Trolls usually pretend to be BLM activists, aiming to divide the democratic party (in this context, being pro-Bernie so that votes are taken away from Hillary). Right trolls imitate Trump supporters, and News Feed handles are “local news aggregators,” typically linking to legitimate news.

For our upcoming analyses, some important variables are:

• author (handle sending the tweet)

• content (text of the tweet)

• language (language of the tweet)

• publish_date (date and time the tweet was sent)

  1. Read in Troll Tweets Dataset

troll_tweets <- read_csv("https://raw.githubusercontent.com/fivethirtyeight/russian-troll-tweets/master/IRAhandle_tweets_12.csv")

2. Basic Data Cleaning and Exploration

1. Remove rows where the tweet was in a language other than English
2. Report the dimensions of the dataset
3. Create two or three basic exploratory plots of the data (ex. plot of the different locations from which tweets were posted, plot of the account category of a tweet)

troll_tweets_clean <- troll_tweets %>% 
  filter(language == "English")

dim(troll_tweets_clean)

## [1] 175966     21

troll_tweets_clean %>% 
  filter(account_type == "Left") %>% 
  ggplot() +
  geom_point(aes(x = following, y = followers, color = updates), alpha = .5)

troll_tweets_clean %>% 
  filter(account_type == "Right") %>% 
  ggplot() +
  geom_point(aes(x = following, y = followers, color = updates), alpha = .5)

troll_tweets_clean %>% 
  filter(account_type == "Russian") %>% 
  ggplot() +
  geom_point(aes(x = following, y = followers, color = updates), alpha = .5)

3. Unnest Tokens: We want each row to represent a word from a tweet, rather than an entire tweet.

troll_tweets_untoken <- troll_tweets_clean %>%
  unnest_tokens(word, content)

troll_tweets_untoken

4. Remove stopwords

troll_tweets_cleaner <- troll_tweets_untoken %>%
  anti_join(stop_words)

troll_tweets_cleaner <- troll_tweets_cleaner %>%
  filter(!word %in% c("http", "https", "t.co", "rt", "amp", 0:9, "a:z"))

5. See how often top words appear

troll_tweets_small <- troll_tweets_cleaner %>%
  count(word) %>% 
  slice_max(order_by = n, n = 50) # 50 most occurring words

# visualize the number of times the 50 top words appear
ggplot(troll_tweets_small, 
       aes(y = fct_reorder(word,n), x = n)) +
  geom_col()

6. Sentiment Analysis

1. Get the sentiments using the “bing” parameter (which classifies words into “positive” or “negative”)
2. Report how many positive and negative words there are in the dataset. Are there more positive or negative words, and why do you think this might be?

get_sentiments("bing")

troll_tweets_sentiment <- troll_tweets_cleaner %>%
  inner_join(sentiments)

troll_tweets_sentiment %>% 
  count(sentiment)

I believe that there are many more negative words because people tend to focus more on negative events, and tweets by bots or trolls tend to only focus on negative events in an attempt to stir up negativity.

7. Using the troll_tweets_small dataset, make a wordcloud:

1. That is sized by the number of times that a word appears in the tweets
2. That is colored by sentiment (positive or negative)

troll_tweets_small %>%
  with(wordcloud(word, n, max.words = 50))

troll_tweets_sentiment %>%
  group_by(word) %>% 
  mutate(n = n()) %>% 
  acast(word ~ sentiment, value.var = "n", fill = 0) %>%
  comparison.cloud(colors = c("red","green"),
                   max.words = 50)

“Undoing” Bias

• The bias starting in the data point is a very concerning point because data collection is so often unfair. However, it seems like a lot of this is due to how data is regulated. From what I understand, all race and gender data is reliant on self-report measures that people are able to opt out of. I also understand that in some situations it is actually illegal to collect race or gender data. Because of this information, I am worried about the future of gender and racial bias. I have always had worries about self-report methods, as they are often unreliable for a multitude of reasons. In certain situations, in order to create fair algorithms, laws on gender and race reporting will have to change for these algorithms to make a true positive difference.